/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/*! \file 
  \brief This file contains the concrete grid
  classes mcrx::dynamic_grid and mcrx::octogrid.

 These classes provide an actual storage for grid cells and override
 the virtual functions get_cell () and the iterator creators begin
 ()/end ().  There is no additional functionality in these classes,
 they only exist because they store their cells in different ways.

 The main point of this whole structure was that since we do grid
 refinement in 2*2*2,we know we will have eight cells and we can avoid
 an extra level of indirection and memory allocation by putting the
 cells directly into the grid object.  Hence, octogrid.  */

// $Id$

// (This file is included automatically by grid.h)


#include "blitz/array.h"

// *** class dynamic_grid ***

/** dynamic_grid is a class containing a dynamically allocated and
    thus arbitrarily sized grid. This is typically the top-level
    grid. */
template <typename cell_data_type,typename sub_grid_type>
class mcrx::dynamic_grid: public grid_base<cell_data_type>  {
public:
  // Redefine to disambiguate the base class types
  typedef grid_base<cell_data_type> T_base;
  typedef typename T_base::T_grid T_grid;
  typedef typename T_base::T_cell T_cell;
  typedef typename T_base::T_cell_tracker T_cell_tracker;
  typedef typename T_base::T_data T_data;

  typedef sub_grid_type T_subgrid;

protected:
  /// The grid cells are stored in a blitz array.
  blitz::Array <grid_cell<cell_data_type> , 3 > cells;

  /// "Virtual constructor" used for creating appropriate sub grids.
  virtual T_grid* create_sub (T_cell *c) const {
    return new sub_grid_type (c);
  };
  /// Placement new versios of the "virtual constructor".
  virtual T_grid* create_sub (T_cell *c, void*& placement) const {
    void* const p = placement;
    reinterpret_cast<sub_grid_type*&> (placement) ++;
    return new (p) sub_grid_type (c);
  };
  void initialize_cells (); ///< Constructor convenience routine.

public:
  dynamic_grid (const vec3d & mi, const vec3d & ma, const vec3i & nn);
  // Creates a dynamic_grid using information from a gridstructure HDU.
  dynamic_grid (CCfits::ExtHDU& file, const vec3d& =vec3d(0,0,0));
  /** Creates an adaptive grid with specified base dimensions and uses
      the grid_factory object to create the grid structure and cell
      data. */
  dynamic_grid (const vec3d & mi, const vec3d & ma, const vec3i & nn,
		const grid_factory<T_data>&); 
  explicit dynamic_grid (T_cell *);
  ~dynamic_grid () {};

  const vec3i& getn () const {
    return cells.shape();};
  
  virtual T_cell *get_cell (const typename T_cell::T_index& indx) {
    return cells.dataZero()+dot (indx, cells.stride());
  }
  virtual const T_cell *get_cell (const typename T_cell::T_index& indx) const {
    return cells.dataZero()+dot (indx, cells.stride());
  }

  // begin () uses the constructor that ensures we step down to a leaf
  // cell if the first cell is refined
  virtual typename T_grid::const_iterator begin () const {
    return typename T_grid::const_iterator (cells.dataFirst ());};
  virtual typename T_grid::iterator begin () {
    return typename T_grid::iterator (cells.dataFirst());};
  virtual typename T_grid::const_local_iterator lbegin () const {
    return typename T_grid::const_local_iterator (cells.dataFirst());};
  virtual typename T_grid::local_iterator lbegin () {
    return typename T_grid::local_iterator (cells.dataFirst());};
  /** Returns an end iterator.  Note that refining cells in the grid
      one level above this (specifically, the cell which follows the
      cell of this grid is in) invalidates this iterator.  */
  virtual typename T_grid::const_iterator end () const {
    if (this->parent_cell()) {
      // if we have a parent cell, the end iterator is in the super grid
      const T_cell* const c = cells.dataFirst() + cells.size ()- 1;
      const typename T_grid::const_iterator e = ++
      (typename T_grid::const_iterator (this, c));
      return e;
    }
    else {
      // if we don't have a parent cell, the end iterator is simply
      // one-past-last
      const T_cell* const c = cells.dataFirst() + cells.size ();
      const typename T_grid::const_iterator e = 
	(typename T_grid::const_iterator (this, c));
      return e;
    }};
  /** Returns an end iterator.  Note that refining cells in the grid
      one level above this (specifically, the cell which follows the
      cell of this grid is in) invalidates this iterator.  */
  virtual typename T_grid::iterator end () {
    if (this->parent_cell()) {
      T_cell* const c = cells.dataFirst() + cells.size ()- 1;
      const typename T_grid::iterator e = ++
	(typename T_grid::iterator (this, c));
      return e;
    }
    else {
      T_cell* const c = cells.dataFirst() + cells.size ();
      const typename T_grid::iterator e = 
	(typename T_grid::iterator (this, c));
      return e;
    }};
  virtual typename T_grid::const_local_iterator lend () const {
    return typename T_grid::const_local_iterator (cells.dataFirst() + 
						  cells.size ());};
  virtual typename T_grid::local_iterator lend () {
    return typename T_grid::local_iterator (cells.dataFirst() + 
					    cells.size ());};
};

// *** class octogrid ***

/** A grid class containing a 2^3 grid. This is used for the subgrids,
    and the fixed size means that we can use statically allocated
    storage for the cells, improving locality of reference and runtime
    efficiency. */
template <typename cell_data_type >
class mcrx::octogrid: public grid_base<cell_data_type>  {
public:
  // Redefine to disambiguate the base class types
  typedef grid_base<cell_data_type> T_base;
  typedef typename T_base::T_grid T_grid;
  typedef typename T_base::T_cell T_cell;
  typedef typename T_base::T_cell_tracker T_cell_tracker;
  typedef typename T_base::T_data T_data;

  typedef octogrid T_subgrid;
protected:

  /** The cells are stored in a tiny vector since we know there are 8.
  
      IMPORTANT: because the tinyvector::operator[] sometimes returns
      COPIES of the cell (the const version does), do NOT use this
      operator to access cells since this will be slow and can even
      lead to incorrect code if we, for example, get a pointer to that
      cell. */
  blitz::TinyVector <T_cell, 8> cells;

  static const vec3i n; ///< This is just so we can return a reference to n.
  
  /// "Virtual constructor" used for creating appropriate sub grids.
  virtual T_grid* create_sub (T_cell *c) const {
    return new octogrid (c);
  };
  /// Placement new version of the virtual constructor.
  virtual T_grid* create_sub (T_cell *c, void*& placement) const {
    void* const p = placement;
    reinterpret_cast<octogrid*&> (placement) ++;
    return new (p) octogrid (c);
  };
  
public:
  // Creates an octogrid as a sub grid in the specified cell.
  explicit octogrid (T_cell*);
  ~octogrid () {};

  virtual const vec3i& getn () const {return n;};
  virtual const T_cell* index (const vec3d& p, bool accept_outside) const;
  virtual T_cell_tracker locate (const vec3d & p, bool accept_outside = false) const;
  
  virtual T_cell *get_cell (const typename T_cell::T_index& indx) {
    return cells.data() +dot (indx, typename T_cell::T_index (4,2,1));
  }
  virtual const T_cell *get_cell (const typename T_cell::T_index& indx) const {
    return cells.data() +dot (indx, typename T_cell::T_index (4,2,1));
  }
  
  virtual typename T_grid::const_iterator begin ()  const {
    return typename T_grid::const_iterator (cells.dataFirst());};
  virtual typename T_grid::iterator begin () {
    return typename T_grid::iterator (cells.dataFirst());};
  virtual typename T_grid::const_local_iterator lbegin ()  const {
    return typename T_grid::const_local_iterator (cells.dataFirst());};
  virtual typename T_grid::local_iterator lbegin () {
    return typename T_grid::local_iterator (cells.dataFirst());};
  /** Returns an end iterator.  Note that refining cells in the grid
      one level above this (specifically, the cell which follows the
      cell of this grid is in) invalidates this iterator.  */
  virtual typename T_grid::const_iterator end () const {
    if (this->parent_cell())
      return ++(typename T_grid::const_iterator (this, cells.data() + 
						 cells.length ()- 1));
    else {
      const T_cell* const c = cells.dataFirst() + cells.length ();
      const typename T_grid::const_iterator e = 
	(typename T_grid::const_iterator (this, c));
      return e;
    };};
  /** Returns an end iterator.  Note that refining cells in the grid
      one level above this (specifically, the cell which follows the
      cell of this grid is in) invalidates this iterator.  */
  virtual typename T_grid::iterator end () {
    if (this->parent_cell())
      return ++(typename T_grid::iterator (this, cells.data() + 
					   cells.length ()- 1));
    else {
      T_cell* const c = cells.dataFirst() + cells.length ();
      const typename T_grid::iterator e = 
	(typename T_grid::iterator (this, c));
      return e;
    };};
  virtual typename T_grid::const_local_iterator lend () const {
    return typename T_grid::const_local_iterator (cells.data() + 
						  cells.length ());};
  virtual typename T_grid::local_iterator lend () {
    return typename T_grid::local_iterator (cells.data() + 
					    cells.length ());};
};

template<typename cell_data_type>
const mcrx::vec3i mcrx::octogrid<cell_data_type>::n (2,2,2);

// *** constructors/destructors ***

/** Creates a dynamic_grid of the specified size. */
template <typename cell_data_type,typename sub_grid_type>
mcrx::dynamic_grid<cell_data_type, sub_grid_type>::
dynamic_grid (const vec3d & mi, const vec3d & ma,
	      const vec3i & nn):
  grid_base<cell_data_type> (mi, ma, nn), cells (nn, blitz::contiguousArray)
{
  assert(cells.isStorageContiguous());
  initialize_cells ();
}

/** Creates a dynamic_grid as a subgrid of another grid. A
    subgrid of this type duplicates the resolution of the parent grid
    (not necessarily so useful, but it's the natural choice). */
template <typename cell_data_type,typename sub_grid_type>
mcrx::dynamic_grid<cell_data_type, sub_grid_type>:: 
dynamic_grid(T_cell *p ):
  grid_base<cell_data_type> (p,p->host_grid ()->getn ()), 
  cells (p->host_grid ()->getn(), blitz::contiguousArray) 
{ 
  assert(cells.isStorageContiguous());
  initialize_cells (); 
}

template <typename cell_data_type,typename sub_grid_type>
void mcrx::dynamic_grid<cell_data_type, sub_grid_type>::initialize_cells ()
{
  typename T_cell::T_index i;
  for (i [0] = 0; i [0] < cells.shape() [0]; ++i [0]) 
    for (i [1] = 0; i [1] < cells.shape() [1]; ++i [1]) 
      for (i [2] = 0; i [2] < cells.shape() [2]; ++i [2]) {
	// This painful expression involving i is because it seems to
	// be the only way I can convert a tiny vector<char> into
	// vec3i...
	T_cell temp(this,i);
	assign_cell (cells (/*cast<int> (i)*/
			    vec3i(i*vec3i(1))), temp);
      }
}


/** Creates an octogrid as a sub grid in a cell. */
template < typename cell_data_type >
mcrx::octogrid<cell_data_type>::
octogrid (mcrx::octogrid<cell_data_type>::T_cell *c ):
  grid_base<cell_data_type> (c, vec3i (2,2,2))
  /*
    For some reason, gcc 4.0 thinks this construct is illegal because
    the grid_cell copy constructor is private.  To me, it's not clear
    the copy constructor is needed anymore in this situation than the
    way it's done in the body.

cells (T_cell (this, typename T_cell::T_index (0,0,0)),
	 T_cell (this, typename T_cell::T_index (0,0,1)),
	 T_cell (this, typename T_cell::T_index (0,1,0)),
	 T_cell (this, typename T_cell::T_index (0,1,1)),
	 T_cell (this, typename T_cell::T_index (1,0,0)),
	 T_cell (this, typename T_cell::T_index (1,0,1)),
	 T_cell (this, typename T_cell::T_index (1,1,0)),
	 T_cell (this, typename T_cell::T_index (1,1,1))) 
  */
{
  {
    T_cell temp (this, typename T_cell::T_index (0,0,0));
    assign_cell(cells[0], temp);
  }
  {
    T_cell temp (this, typename T_cell::T_index (0,0,1));
    assign_cell(cells[1], temp);
  }
  {
    T_cell temp (this, typename T_cell::T_index (0,1,0));
    assign_cell(cells[2], temp);
  }
  {
    T_cell temp (this, typename T_cell::T_index (0,1,1));
    assign_cell(cells[3], temp);
  }
  {
    T_cell temp (this, typename T_cell::T_index (1,0,0));
    assign_cell(cells[4], temp);
  }
  {
    T_cell temp (this, typename T_cell::T_index (1,0,1));
    assign_cell(cells[5], temp);
  }
  {
    T_cell temp (this, typename T_cell::T_index (1,1,0));
    assign_cell(cells[6], temp);
  }
  {
    T_cell temp (this, typename T_cell::T_index (1,1,1));
    assign_cell(cells[7], temp);
  }
}


/** Returns the index of the cell containing a specified
    position. This is a redefinition of the grid_base version of
    index().  Since we know that we are in a 2^3 grid, we can be much
    more efficient. */
template < typename cell_data_type >
inline const typename mcrx::octogrid<cell_data_type>::T_cell*
mcrx::octogrid<cell_data_type>::index (const vec3d& p,
				       bool accept_outside) const
{
  //  calculate "floating point index" (it can be between -1 and 0 so
  //  we have to use floating points to detect this since it would
  //  round to 0.
  const vec3d iii=(p- this->getmin() )*this->getinvsize ();

  // Calculate which cell corresponds to this, assuming that we accept
  // outside of the grid
  int i = ((iii[0] > 1) << 2) | ((iii[1] > 1) << 1) | ((iii[2] > 1));
  assert (i < 8);
  assert (i >= 0);
  
  // Do we have to do the outside check? Why is the parent cell check
  // here?  Is this because if we've been called from higher up they
  // have already decided that the point should be in this grid?
  if (!accept_outside && !this->parent_cell () &&
      // Yes, the check is required
      (any (iii < 0) || any (iii > 2)))
    // We are outside of the grid
    return 0; 
  else
    // No, we can just return
    return cells.data() +i;
}


/** Find the grid cell that contains position p.  If this grid is not
 in a top-level grid, we accept positions outside of the grid itself,
 and just use the closest coordinates. This is used by neighbor. */
template < typename cell_data_type >
inline typename mcrx::octogrid<cell_data_type>::T_cell_tracker
mcrx::octogrid<cell_data_type>::locate(const vec3d& p,
				       bool accept_outside) const
{
  // Now we can make a nonvirtual call here
  const T_cell* const c = octogrid<cell_data_type>::index (p, accept_outside);
  
  if(c && !c->is_leaf())
    // That cell exists but is subdivided, recurse in the locate call
    // (accept_outside is not necessary in this case because it only
    // affects the top-level grid)
    
    // The reinterpret cast is a hack to avoid calling locate
    // virtually -- we KNOW that the sub grid of an octogrid is itself
    // an octogrid. Ugly, but you do what you have to...
    return reinterpret_cast<const octogrid<cell_data_type>*>
      (c->sub_grid()) ->octogrid<cell_data_type>::locate (p);
  else
    // cell is not subdivided, just return it
    return const_cast<T_cell*>(c);
}
